---
title: "Alternative coding of democracy for 'Can courts in non-democracies deter election fraud? De jure judicial independence, political competition, and election integrity'"
output: html_notebook
author: Cole J. Harvey
date: 12/28/2021
---

```{r setup}
library(tidyverse)
library(lme4)
library(lmtest)
library(ggplot2)
library(ggeffects)
library(interplot)
library(interactions)
library(interplot.medline)
library(interflex)
library(stargazer)
library(wfe)
library(plm)
library(ebal)
library(cowplot)
library(sjPlot)
library(fastDummies)



#####Loading dataset####
vdem.nodems <- read.csv("harvey_courts_manipulation_appendix_data_altdem.csv")
vdem.nodems <- vdem.nodems %>% mutate(COWyear = paste(COWcode, year, sep = "_"))


vdem.nodems.elections <- subset(vdem.nodems, is.na(vdem.nodems$v2elirreg.inv)==FALSE)


###Coding presidential elections as majoritarian
vdem.nodems.elections <- vdem.nodems.elections %>% mutate(v2elparlel = ifelse(elexec == 1, 0, v2elparlel))

##Dummy variables for electoral system
vdem.nodems.elections <- fastDummies::dummy_cols(vdem.nodems.elections, select_columns = "v2elparlel")

vdem.nodems.elections$COW.factor <- factor(vdem.nodems.elections$COWcode)
vdem.nodems.elections$year.factor <- factor(vdem.nodems.elections$year)

vdem.nodems.elections <- vdem.nodems.elections %>% filter(parcomp.lag >= 0)
vdem.nodems.elections <- vdem.nodems.elections %>% filter(is.na(reform_positive.lag) == F)

##Lagging election variables

vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2elirreg.inv.lag = dplyr::lag(v2elirreg.inv))
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2elintim.inv.lag = dplyr::lag(v2elintim.inv))   
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2eldommon.lag = dplyr::lag(v2eldommon)) 
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(elexec.lag = dplyr::lag(elexec))   


##Leading election variables

vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2elirreg.inv.lead = dplyr::lead(v2elirreg.inv))
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2elintim.inv.lead = dplyr::lead(v2elintim.inv))   
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2eldommon.lead = dplyr::lead(v2eldommon)) 
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(elexec.lead = dplyr::lead(elexec))   
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(opposition.oversight.lead = dplyr::lead(v2lgoppart))   
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(gdp.lead = dplyr::lead(e_migdppcln))
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2elintmon.lead = dplyr::lead(v2elintmon))  
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2elparlel_1.lead = dplyr::lead(v2elparlel_1))  
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(COWcode) %>% mutate(v2elparlel_2.lead = dplyr::lead(v2elparlel_2))  


##Needed variables from vdem

vdem <- vdemdata::vdem
vdem <- vdem %>% dplyr::select(v2x_regime, v2x_clphy, e_migdpgro, v2x_cspart, v2elvotlrg, country_id, year)
vdem <- vdem %>% dplyr::filter(year >= 1944)
vdem <- vdem %>% mutate(COWcode = countrycode::countrycode(sourcevar = country_id, origin = "vdem", destination = "cown"))

vdem <- vdem %>% mutate(COWyear = paste(COWcode, year, sep = "_"))
vdem.nodems.elections <- vdem.nodems.elections %>% left_join(vdem, by = "COWyear")

vdem.nodems.elections <- vdem.nodems.elections %>% mutate(winner.margin = ifelse(elexec == 1, v2elvotlrg/100, v2ellovtlg/100))
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(country_id) %>% mutate(cspart.1lag = lag(v2x_cspart))
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(country_id) %>% mutate(gdpgro.1lag = lag(e_migdpgro))
vdem.nodems.elections <- vdem.nodems.elections %>% group_by(country_id) %>% mutate(physinteg.1lag = lag(v2x_clphy))

vdem.nodems.elections <- fastDummies::dummy_cols(vdem.nodems.elections, select_columns = "v2x_regime")

vdem.nodems.elections <- vdem.nodems.elections %>% filter(v2x_regime_3 == 0)


vdem.nodems.elections <- vdem.nodems.elections %>% dplyr::select(-year.y)
vdem.nodems.elections <- vdem.nodems.elections %>% rename(year = year.x)

```

##Main models


Elcomp

```{r ebal elcomp, include=FALSE}



myvars <- c("COW.factor", "year", "transitional", "v2elirreg.inv", "v2elintim.inv", "lc.ind.2lag",  "hc.ind.2lag", "oppaut.2lag", "elexec", "education.2lag",  "v2elparlel_1", "v2elparlel_2", "lji.2lag", "v2elparlel_3", "purge.lag", "pack.lag",
            "elcompindex.m.lag", "elcompindex.m.2lag", "urban.2lag",  "v2elintmon", 
            "country_name", "reform_positive.lag", "loggpdpc.2lag", "reform_negative.lag",
            "democracy.duration.2lag", "exec.respectcon.2lag", "leg.constraints.2lag", 
            "altinfo.2lag", "loggpdpc.lag") 

dataset.matching <- vdem.nodems.elections[myvars]
dataset.matching.complete <- na.omit(dataset.matching)


myvars <- c("democracy.duration.2lag", "oppaut.2lag",  "elcompindex.m.2lag", "lji.2lag",
             "loggpdpc.2lag", "urban.2lag",
            "hc.ind.2lag",  "lc.ind.2lag", "transitional", "exec.respectcon.2lag", "altinfo.2lag", "education.2lag",
            "leg.constraints.2lag") 
ebal.covariates <- dataset.matching.complete[myvars]

ebal.test <- ebalance(Treatment=dataset.matching.complete$reform_positive.lag, X = ebal.covariates)



# means in treatment group data
t.means <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==1,],2,mean)
# means in reweighted control group data
c.means.bal <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==0,],2,weighted.mean,w=ebal.test$w)
# means in raw data control group data
c.means <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==0,],2,mean)



dataset.matching.complete.controls <- subset(dataset.matching.complete, dataset.matching.complete$reform_positive.lag == 0)
dataset.matching.complete.controls <- data.frame(cbind(dataset.matching.complete.controls, ebal.test$w))


dataset.matching.complete.treat <- subset(dataset.matching.complete, dataset.matching.complete$reform_positive.lag == 1)
dataset.matching.complete.treat$ebal.test.w <- 1

dataset.matching.complete.w <- bind_rows(dataset.matching.complete.treat, dataset.matching.complete.controls)


###regression
elirreg.posreform.elcomp.base <- lm(v2elirreg.inv~reform_positive.lag +
                                        
                                        elcompindex.m.lag  + elexec +  v2elparlel_1 + v2elparlel_2 + 
                                        loggpdpc.lag + 
                                            v2elintmon +
                                            reform_negative.lag + purge.lag + pack.lag+
                                      COW.factor, 
                                    weights=ebal.test.w,  
                                    data = dataset.matching.complete.w)
#summary(elirreg.posreform.elcomp.base)


mm.elirreg.posreform.elcomp.all <- lm(v2elirreg.inv~reform_positive.lag
                                        + 
                                          elcompindex.m.lag +
                                      reform_positive.lag*elcompindex.m.lag
                                      + elexec + v2elparlel_1 + v2elparlel_2 + 
                                          loggpdpc.lag + 
                                              v2elintmon +
                                          reform_negative.lag+ purge.lag + pack.lag+
                                     COW.factor,  
                                            weights=ebal.test.w,
                                            data = dataset.matching.complete.w)
#summary(mm.elirreg.posreform.elcomp.all)

lrtest(elirreg.posreform.elcomp.base, mm.elirreg.posreform.elcomp.all)

```


```{r plot ebal elcomp}
p2.elcomp.ch <- interplot(mm.elirreg.posreform.elcomp.all, var1="reform_positive.lag", var2="elcompindex.m.lag",
                hist=TRUE, adjCI = TRUE) + theme_bw() + labs(x="Political openness \n(1-year lag)", y="Marginal effect") +
  geom_hline(yintercept=0, linetype=2)  
p2.elcomp.ch
```

Legislative oversight

```{r include=FALSE}
### Legislative opposition oversight
myvars <- c("COW.factor", "year.factor", "transitional", "v2elirreg.inv", "v2elintim.inv", "lc.ind.2lag", 
            "hc.ind.2lag", "oppaut.2lag", "elexec", "education.2lag",  "lji.2lag", "pack.lag", "purge.lag",
            "opposition.oversight.lag", "opp.oversight.2lag", "urban.2lag", "v2elmulpar", "v2elintmon", 
            "country_name", "reform_positive.lag", "loggpdpc.2lag", "reform_negative.lag",
            "democracy.duration.2lag", "exec.respectcon.2lag", "leg.constraints.2lag", 
            "altinfo.2lag", "loggpdpc.lag", "years.since.election",  "v2elparlel_1", "v2elparlel_2",  "v2elparlel_3")  # 

dataset.matching <- vdem.nodems.elections[myvars]
dataset.matching.complete <- na.omit(dataset.matching)


myvars <- c("democracy.duration.2lag", "oppaut.2lag",  "opp.oversight.2lag", "lji.2lag",
            "loggpdpc.2lag", "urban.2lag",
            "hc.ind.2lag",  "lc.ind.2lag", "transitional", "exec.respectcon.2lag", "altinfo.2lag", "education.2lag",
            "leg.constraints.2lag") 
ebal.covariates <- dataset.matching.complete[myvars]

ebal.test <- ebalance(Treatment=dataset.matching.complete$reform_positive.lag, X = ebal.covariates)




# means in treatment group data
t.means <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==1,],2,mean)
# means in reweighted control group data
c.means.bal <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==0,],2,weighted.mean,w=ebal.test$w)
# means in raw data control group data
c.means <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==0,],2,mean)


dataset.matching.complete.controls <- subset(dataset.matching.complete, dataset.matching.complete$reform_positive.lag == 0)
dataset.matching.complete.controls <- data.frame(cbind(dataset.matching.complete.controls, ebal.test$w))


dataset.matching.complete.treat <- subset(dataset.matching.complete, dataset.matching.complete$reform_positive.lag == 1)
dataset.matching.complete.treat$ebal.test.w <- 1

dataset.matching.complete.w <- bind_rows(dataset.matching.complete.treat, dataset.matching.complete.controls)


###regression
elirreg.posreform.oversight.base <- lm(v2elirreg.inv~reform_positive.lag +
                                               opposition.oversight.lag + elexec + 
                                           loggpdpc.lag + 
                                               v2elintmon + 
                                            v2elparlel_1 + v2elparlel_2 +  reform_negative.lag + purge.lag + pack.lag+
                                         COW.factor 
                                             , weights=ebal.test.w,
                                             data = dataset.matching.complete.w)
#summary(elirreg.posreform.oversight.base)


mm.elirreg.posreform.oversight <- lm(v2elirreg.inv~reform_positive.lag +
                                             opposition.oversight.lag + reform_positive.lag*opposition.oversight.lag   + elexec + 
                                         loggpdpc.lag + 
                                              v2elintmon +
                                              v2elparlel_1 + v2elparlel_2   + reform_negative.lag + purge.lag + pack.lag+
                                      COW.factor, 
                                           weights=ebal.test.w,
                                           data = dataset.matching.complete.w)
#summary(mm.elirreg.posreform.oversight)

```


```{r}
p2.legov.ch <- interplot(mm.elirreg.posreform.oversight, var1="reform_positive.lag", var2="opposition.oversight.lag",
                hist=TRUE, adjCI = TRUE) + theme_bw() + labs(x="Opposition oversight \n(1-year lag)", y="Marginal effect") +
  geom_hline(yintercept=0, linetype=2)   
p2.legov.ch
```



Polcon


```{r ebal polcon, include=FALSE}
vdem.nodems.elections <- vdem.nodems.elections %>% mutate(scale.polcon.lag = log(h_polcon3.lag))
myvars <- c("COW.factor", "year.factor", "transitional", "v2elirreg.inv", "v2elintim.inv", "lc.ind.2lag", "lji.2lag", "hc.ind.2lag", "oppaut.2lag", "elexec", "education.2lag",  "v2elparlel_1", "v2elparlel_2",                        "h_polcon3.lag", "h_polcon3.2lag",
             "urban.2lag",  "v2elintmon", "pack.lag", "purge.lag",
            "country_name", "reform_positive.lag", "loggpdpc.2lag", "reform_negative.lag",
            "democracy.duration.2lag", "exec.respectcon.2lag", "leg.constraints.2lag", 
            "altinfo.2lag", "loggpdpc.lag") 

dataset.matching <- vdem.nodems.elections[myvars]
dataset.matching.complete <- na.omit(dataset.matching)


myvars <- c("democracy.duration.2lag", "oppaut.2lag",  "h_polcon3.2lag", "lji.2lag",
             "loggpdpc.2lag", "urban.2lag",
            "hc.ind.2lag",  "lc.ind.2lag", "transitional", "exec.respectcon.2lag", "altinfo.2lag", "education.2lag",
            "leg.constraints.2lag")  
ebal.covariates <- dataset.matching.complete[myvars]

ebal.test <- ebalance(Treatment=dataset.matching.complete$reform_positive.lag, X = ebal.covariates)



# means in treatment group data
t.means <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==1,],2,mean)
# means in reweighted control group data
c.means.bal <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==0,],2,weighted.mean,w=ebal.test$w)
# means in raw data control group data
c.means <- apply(ebal.covariates[dataset.matching.complete$reform_positive.lag==0,],2,mean)



dataset.matching.complete.controls <- subset(dataset.matching.complete, dataset.matching.complete$reform_positive.lag == 0)
dataset.matching.complete.controls <- data.frame(cbind(dataset.matching.complete.controls, ebal.test$w))


dataset.matching.complete.treat <- subset(dataset.matching.complete, dataset.matching.complete$reform_positive.lag == 1)
dataset.matching.complete.treat$ebal.test.w <- 1

dataset.matching.complete.w <- bind_rows(dataset.matching.complete.treat, dataset.matching.complete.controls)


###regression
elirreg.posreform.polcon.base <- lm(v2elirreg.inv~reform_positive.lag +
                                        
                                        h_polcon3.lag  + elexec +  v2elparlel_1 + v2elparlel_2  +
                                                v2elintmon +
                                           reform_negative.lag + loggpdpc.lag + purge.lag + pack.lag+
                                      COW.factor, 
                                    weights=ebal.test.w,  
                                    data = dataset.matching.complete.w)
#summary(elirreg.posreform.polity.base)


mm.elirreg.posreform.polcon.all <- lm(v2elirreg.inv~reform_positive.lag
                                        + 
                                          h_polcon3.lag +
                                      reform_positive.lag*h_polcon3.lag
                                      + elexec + v2elparlel_1 + v2elparlel_2  +
                                                  v2elintmon  +
                                               reform_negative.lag + loggpdpc.lag + purge.lag + pack.lag+
                                     COW.factor,  
                                            weights=ebal.test.w,
                                            data = dataset.matching.complete.w)
#summary(mm.elirreg.posreform.polcon.all)


```

```{r}
p4.polcon.ch <- interplot(mm.elirreg.posreform.polcon.all, var1="reform_positive.lag", var2="h_polcon3.lag",
                hist=TRUE, adjCI = TRUE) + theme_bw() + labs(x="Political constraints (1-year lag)", y="Marginal effect") +
  geom_hline(yintercept=0, linetype=2)   
p4.polcon.ch
```

## Producing Figure B.5

```{r combining me plots}
p.combined.ebal <- cowplot::plot_grid(p2.elcomp.ch, p2.legov.ch, p4.polcon.ch, align = "hv", nrow = 2, label_size = 12, hjust = 0)
p.combined.ebal

```


## Producing Table B.5

```{r}
#Entropy balanced
stargazer(mm.elirreg.posreform.elcomp.all, mm.elirreg.posreform.oversight, mm.elirreg.posreform.polcon.all,
          digits = 2, omit = c("COW.factor", "COWcode", "year.factor", "year"),
          type="text")


```